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Abstract 



Passive molecular movements in cells are often claimed to exhibit "anoma- 
lous subdiffusion" , where the mean-squared displacement (MSD) scales as 
a power law of time with exponent a < I. Diffusion hindrance by ob- 
stacles is often invoked to explain these observations. In many studies of 
hindered diffusion, the estimated values of a vary strongly. This led to the 
hypothesis that a depends on obstacle density. This is however at odds 
with the theoretical support for hindered diffusion among randomly located 
immobile obstacles, which predicts that true subdiffusion occurs only in the 
vicinity of the threshold for the percolation of obstacles, and that a takes 
a unique, universal value. Here, we present refined simulations of hindered 
diffusion with biologically realistic parameters and bring forth four main 
contributions, (i) We confirm that diffusion hindered by randomly located 
immobile obstacles does not exhibit variations of a and (ii) that the MSD 
in fact never scales like a power law of time, (in) In contrast to diffusing 
obstacles, obstacles fluctuating around equilibrium positions preserve and 
even emphasize anomalous regimes, (iv) Hindered diffusion is not equiva- 
lent to anomalous diffusion due to random traps with heavy-tailed trap time 
distribution. These results shed new light on the existing literature about 
subdiffusion. 

Key words: Anomalous Diffusion; Obstacle Hindered Diffusion; Macro- 
molecular Crowding; Fluctuating Obstacles; Monte-Carlo Simulations; Per- 
colation Theory 



Introduction 



A hallmark of Brownian motion is the linear relation between the average of 
the squared distance travelled by the molecule (mean squared displacement) 
and time: oc t. However, experimental measurements of molecular 

diffusion in living cells have consistently reported nonlinear relations in al- 
most all cell compartments, either in procaryotes or eucaryotes. Most often, 
these nonlinear variations are claimed to be power laws {{R^{t)) (X r with 
a 7^ 1), i.e. straight lines once represented on log-log plots. 

Super-diffusive motion (2 > a > 1) is relatively well understood, being 
usually due to active transport mediated by molecular motors on cytoskele- 
ton elements (jll), but subdiffusive transport {a < 1) less so. Fundamen- 
tally, it is expected to arise when a hierarchy of slow processes hinders mo- 
tion. Two classes of mechanisms yield such a situation: (a) a hierarchy of 
time scales exists due to transient capture of the diffusive molecule by traps 
(e.g. weak binding) with broadly-distributed trapping times (this is the so- 
called continuous-time random walks, or CTRW, model), or (b) it can result 
from the heterogeneous, irregular geometry of the space domain in which 
molecules move. In measurements in cells, it is usually unknown which one 
of these mechanisms is the cause of the observed anomalous behavior. This 
is a first source of ambiguity in the interpretation of experimental results. 

Moreover, experimental works report strikingly different values of a. 
Continuously- varying exponents are consistent with the CTRW model, where 
the trap time distribution is power-law distributed p{t) oc t~/^ with 1 < /3 < 
2: in this case, diffusion will be anomalous, with exponent a = (3 — 1 
However, it is not clear how such a distribution could arise in vivo. More- 
over, several experiments have reported anomalous diffusion for cases where 
trapping can be excluded, i.e. for purely obstacle-hindered diffusion (macro- 
molecular crowding) (Q,!!)- So, even though trapping can play a role in vivo, 
hindered diffusion has to be considered. 

Theoretically, hindered diffusion, i.e. diffusion in the presence of volume- 
excluding immobile obstacles, corresponds to the movement of a random 
walker on a percolating cluster (the set left available by the obstacles). At 
an obstacle density at finite distance above the percolation threshold, the 
macroscopic percolating cluster is fractal only over short distances and ap- 
pears uniform at larger scales (j6i). One expects then three diffusion regimes: 
Brownian (a = 1) at short times (before the walker encounters obstacles), 
followed by a crossover to an anomalous regime with exponent q < 1 at 
intermediate times and, finally, a second crossover back to normal diffusion 
(a = 1) 0), 

which is the asymptotic regime. When the density of obstacles 
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changes, the two crossover times are predicted to vary: in particular the 
latter (from anomalous back to Brownian) is expected to diverge when the 
percolation threshold is approached (0). But a major prediction of the the- 
ory is that the value of the exponent a is universal, i.e. it does not depend 
on obstacle density, but essentially only on the (embedding) space dimen- 
sion, taking respectively values a ~ 0.686 and ~ 0.526 in two and three 
dimensions ((3). 

Inside bacteria, experimental measurements of a for RNA particles in 
the cytoplasm gave estimates in the [0.7 — 0.8] range (0). In the cytoplasm 
and nucleus of mammalian cells, diffusion of artificial markers yielded values 
ranging from [0.5 - 0.6] ft, to [0.6 - 0.8] ji). Such a large variability was 
also reported for GFP in the cytoplasm of mammalian cells (0). Reports 
of molecule diffusion on the plasma membrane of mammal cells are even less 
coherent. The diffusion of small fluorescent probes led to q ~ 0.74 (i.e. not 
far from the value expected at the percolation threshold for 2d) (11). But 
larger lipid-attached particles yielded a value smaller than this theoretical 
limit (a w 0.53) (fl^). and studies using fluorescent-tagged receptors have 
reported huge variability across cells (a G [0.2-1.0] HH)). Finally, in the 
nucleus of mammal cells, the diffusion of Cajal bodies was found anomalous, 
with a ~ 0.67 (llSl ). But the diffusion of GFP or related proteins in nuclei 
also demonstrated large variability (a G [0.65 — 0.9] (fiol )). 

Such wildly varying estimates leave a confusing situation. To be true, 
most of the experimental curves in the above reports fail to display clearly 
identifiable power-law regimes, e.g. clean straight lines in log- log scales (12i). 
In many cases, the anomalous regime is too short, i.e. lasts hardly more than 
a decade, to allow for the unambiguous identification of a true power-law 
regime and/or the number of points used to fit the data is too small (1^, 15). 

One would expect that the interpretation of experimental measurements 
in vitro should be easier. Diffusion in entangled meshes of semiflexible poly- 
mers can yield anomalous diffusion with exponent = 3/4 for particles 
larger than the mesh size (|l6|). This regime has been observed during the 
diffusion of large particles in actin solutions ([17|) and in the cytosplasm of 
eucaryotic cells (see e.g. (fisl)). The mesh size of the actin network in vivo is 
however of the order of 1 ^m ([itI ) so that this anomalous regime is expected 
for molecular sizes at least several hundreds of nm, circa 2 orders of magni- 
tude larger than most proteins. Diffusion of smaller molecules (< 5 — 6 nm) 
in polymer solutions is thus generally expected to be Brownian (a = 1) but 



with a reduced apparent diffusion coefficient ( 19 ]^ 
Nevertheless, several experimental studies ([91, 



2ol . 21 ) have reported 



anoma- 



lous diffusion of small artificial molecules or globular proteins in crowded 
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dextran solutions, claiming the observation of anomalous exponents with val- 
ues decreasing from 1 to ~ 0.70 when molecular crowding increases (0, [iol ). 
However, experimental error on the determination of a seems very large. 
Repeated measurements of a on the same experimental sample yield values 
that can vary from 0.7 to 0.9 across experiments (j2ll ). These results are fur- 



ther challenged by several contradictory experimental studies (|l9l . |22| . |23| ). 
that found only Brownian diffusion in similar or identical conditions. Taken 
together, these studies attest that even in in vitro and carefully controlled 
conditions, the processing of data is a step not free from various pitfalls, 
and the large variability of the reported results calls for caution in assessing 
their worth. 

More intriguingly, even computer simulations of hindered diffusion have 
claimed crowding-dependent values of a (|24l . l25l). This notion was first sug- 
gested more than 25 years ago in Argyrakis and Kopelman (|2^). where the 
diffusion of a random walker on a percolation cluster above the threshold 
was studied by computer simulations. However, the power law regime at 
intermediate times claimed in this paper is poorly convincing, especially 
because the simulations barely cover two time decades. Similarly, in pub- 
lished simulation results, the existence of a true power-law regime in the 
(log(i?2(t)) vs logt) curves is questionable (at least far from the thresh- 
old). First of all, the size of the spatial domain used in the simulations is 
very often small compared to the maximal average distance travelled by the 
molecules (|24l . |25| . 1271 ) so that finite size effects are to be expected. More- 
over, simulation results often suffer from the same weaknesses as experimen- 
tal ones, including small number of data points ([^21)_^or short anomalous 



regimes with disputable evidence of a power law ([281, l29l ) . 

Finally, it has been suggested that the anomalous regime observed may 
be crucially dependent on the dynamics (or lack of) of obstacles themselves. 
Indeed, when the obstacles are allowed to diffuse even at a much slower rate 
than the tracked molecule, the transient anomalous regime is predicted to 
disappear rapidly (271. l30l). In fact, a recent theoretical study (31') predicts 
that when mobile obstacles and walkers have the same diffusion coefficient, 
the movement would be Brownian diffusion at long times and anomalous at 
short ones. But this predicted short-time anomalous diffusion regime is not 
a power law. 

More recent works have moreover suggested an interesting possibility 
in the case where obstacle diffusion is small but non zero. In this case, one 
can envision the mobile obstacles as dynamic cages, transiently trapping the 
diffusing molecules. The lifetimes of these cages may be broadly distributed 
at high obstacle density, thus potentially leading to a situation similar to the 
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trapping-mediated anomalous diffusion (CTRW) evoked above |32). Sucli a 
situation would offer a plausible explanation for the reported variability of 
the anomalous exponent but remains to be tested by realistic simulations. 

In the present work, our objective is to examine the issues raised above. 
Using refined simulations devoid of the shortcomings evoked above, we make 
four contributions : 

• We confirm that hindered diffusion due to immobile obstacles does not 
give rise to anomalous regimes with crowding-dependent exponents. 
The anomalous exponent is unique, and obstacle density only modifies 
the crossover times between Brownian and anomalous regimes. 

• For biologically realistic parameters, our simulations do not show ev- 
idence for even a transient power law behavior during the transient 
anomalous regime: the dynamics in this regime are completely domi- 
nated by crossovers. 

• The effect of obstacle movements mainly depends on the type of move- 



ment considered. In agreement with previous works (|27l. I3d). our sim- 
ulations show that diffusing obstacles (undergoing Brownian motion) 
have a strong tendency to suppress anomalous regimes. But if the 
obstacles jitter around equilibrium average positions, the anomalous 
regime can persist at short times even with strong obstacle fluctua- 
tions. 

• Finally, we show that, whatever the type of obstacle considered, one 
can measure a distribution of "effective trapping times" which slightly 
broadens with increasing obstacle density, but without exhibiting the 
heavy-tails characteristic of the CTRW model. 



Methods 

Simulations with immobile obstacles 

We simulated the diffusion of tracer molecules using lattice-free domains 
with periodic boundary conditions (this reduces finite-size effects). Each run 
was initiated by positioning at random (with uniform probability) 2d (3d) 
obstacles in a 2d (3d) continuous domain of overall size x Ly{xLz), until 
the volume fraction occupied by the obstacles equals the preset excluded 
volume fraction 9. A set of A^w non- interacting tracer molecules (walkers) 
was then positioned in the space domain at random locations but respecting 
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excluded volume w.r.t. obstacles: for each tracer molecule, a new position 
Xyj is chosen at random (with uniform distribution) inside the simulation 
domain, until the tracer molecule docs not overlap with any obstacle at x^. 

At each time step, the simulation then proceeds by moving each tracer 
molecule (walker) independently of each other. We modeled diffusion using 
values for the time step 5t and the displacement per time step 5r that do 
not depend on the coefficient diffusion. Between t and t + 6t, each molecule 
has a probability Pmove to move to a randomly-chosen position located at 
distance Sr from its position at t. The displacement probability is given 
by the diffusion coefficient D : Pmove = 2dD6t/5r'^ . The advantage of 
this algorithm is that, choosing a sufficiently small value for 5r (namely 
6r < 2r^ + 2robs) ensures the excluded volume condition for all values of 
D and 5t. In our simulations, typical values ra.ng6d. from -Fmove — 

0.35 

to 0.95. Excluded volume is then modeled by adding the restriction that 
displacement attempts are rejected when the diffusing molecule, at the target 
site, overlaps with an obstacle. 

The squared displacement R^{t) of each tracer molecule was monitored 
taking into account periodic boundary conditions. Unless otherwise indi- 
cated, R^{t) was averaged over = 10^ molecule trajectories per spatial 
obstacle configuration and 10 obstacle configurations. 

Simulations with mobile obstacles 

To model the diffusive movements of the obstacles, the position of each 

obstacle at time t + 5t was updated according to robs(i + ^t) = i'obs(i) + 
(cr^r) where (crsr) is a d-dimensional random vector of which each 
component is drawn independently from a normal distribution of mean 
and standard deviation asr- This setting results in a diffusive movement 
with diffusion coefficient Dobs = <^sr/ i^St). 

Likewise, to simulate fluctuating obstacles, the position of each obstacle 
at time t + 6t was set according to robs(i + ^t) = robs(^o) + {asr) where 
robs(io) is the equilibrium position of the obstacle . The quantity asr thus 
represents the amplitude of the obstacle fluctuations around their equilib- 
rium position, i.e. their mobility. Note that in both cases obstacles do not 
experience excluded volume interactions with respect to each other and can 
thus interpenetrate. 
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Effective and approximate exponents 

In many cases, a simple log-log plot of the MSD as a function of 

time is not accurate enough to validate or invalidate the existence of a power 
law regime. To improve the accuracy of the measurements, we computed 
the evolution with time of an effective exponent, a^sit), the local slope of 
the vs t relation in log-log scale: a^sit) = dlog /dlogt. The 

presence of a power-law regime then translates to an easily recognizable 
plateau in the acfr(^) vs t curve. 

Even though most of our results show that the vs t relation does 

not evidence a power law regime (ie a plateau in the aeff(i) vs t curve), we 
define an approximate exponent, a as the minimal value of the aes{t) vs t 
curve, to allow for comparisons with published results. 

Trapping time distribution 

Mobile obstacles at high density are expected to form "dynamic cages" for 
the diffusing tracer molecule, that transiently form and disappear in time. 
Because of their stochastic and transient nature, they are difficult to define 
strictly. Here we consider that a tracer molecule is trapped by an obstacle 
whenever the tracer molecule remains within a distance rtrap (the trapping 
distance) of the obstacle. We then monitor the number of consecutive time 
steps spent by the tracer molecule within rtrap of every obstacle. These 
trapping times r for each encounter of each obstacle and each molecule were 
then pooled to yield the trapping time distribution. 

Results 

Hindered diffusion by immobile obstacles 

Our first goal here is to question the possibility that, for biologically relevant 
parameters, diffusion among immobile obstacles can give rise to anomalous 
scaling regimes with exponents that vary with the density of obstacles. 

To avoid the issues encountered by previous simulation studies of hin- 
dered diffusion, we took care to monitor the mean-square displacement over 
large time scales (at least 6 decades) with a good sampling (> 700 data 
points per curve). The size of the spatial domain was set to Lx = Ly = 
5 X 10^ X Sr in 2d simulations and = Ly = = 10^ x Sr in 3d ones. 
Using these settings, the maximal mean travelled distance varied from 8.4 
to 16.8 % of the length of the spatial domain in 2d, hence excluding finite 
size effects. 
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The typical values of the molecule and obstacle sizes and mobility were 
chosen so as to be representative of the size encountered in a typical bacterial 
cell. We considered an "average" protein in the cytoplasm of E. coli with ra- 
dius = 2.0 nm and molecular weight 40 kDa (see e.g. TableSl in McGuf- 
fee and Elcock (jsi)). The diffusion coefficient of the 28 kDa GFP in E. coli 
cytoplasm ranges from 2 to 10 /im^/s while that of the 72 kDa GFP-MBP 
fusion protein is around 2.5 /xm^/s (j34l .l35l). We accordingly set the diffusion 
coefficient of our 40 kDa "average" protein to within D £ [4.0 — 8.0] ^m^/s. 
We considered large multimolecular obstacles, such as ribosomes, with ra- 
dius Tobs = 5 nm. In 3d simulations however, computational load demanded 
larger 6r values and we had to increase robs to 10.0 nm. 

Figure [1] shows the evolution with time of the MSD (Fig. [T]A.) as a log- log 
plot. Each curve in the figure corresponds to a different obstacle density. 
The top-most curve is for unobstructed diffusion, while obstacle density 
increases from top to bottom, until the lowest curve where 50% of space is 
occupied by obstacles (thus the excluded fraction 9 = 0.50). 

Clearly, without obstruction, a single diffusion regime is found, with 
Brownian diffusion (a = 1). When the excluded fraction is very high (e.g., 
9 = 0.50, well above the percolation threshold located around 9 = 0.44, 
see below), the curves evidence saturation of the MSD, as space is parti- 
tioned into disconnected clusters of available sites which permanently trap 
the tracer molecules. For intermediate obstacle fractions {9 < 0.4), three 
regimes can be distinguished, including two Brownian regimes: one at short 
times {t < 10~^ ms) and another one at long times (t > 0.5 ms). Note that, 
within the time scale of figure [H the late Brownian regime is clearly reached 
only for the smallest obstacle densities. For larger obstacle densities, the 
curves are however clearly converging to it. 

As is often the case, it is very delicate to tell from the (log MSD vs log t) 
curves what happens in the intermediate-time regime. This is better seen 
from the time evolution of the local exponent aeff(t) = d log (-R^(t)) /dlogt 
(see Methods), which is shown in the bottom panel of Fig. [T]A.. Without ob- 
stacles, the local slope is 1.0 at all times. With increasing excluded fraction 
9, Oesit) first bends downward then converges back to 1.0 at long times. For 
the largest excluded fractions, the local slope shows a monotonous tendency 
to decay at long times, corresponding to the saturation seen in the MSD 
vs time curves (Fig. [T]A, Top). The visual change in aspect of the local 
slope curves (upward vs downward curvature at long times) can be used to 
locate the percolation threshold (0). This allows us to estimate the criti- 
cal excluded fraction 9c = 0.441(2) (where the digit in parenthesis indicates 
incertitude). The long-time value of aeff(i) for 9 = 9c yields an estimate for 
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the anomalous exponent a (at the threshold) : a = aefr(i — > oo,^ — > ^c) = 
0.68(4). This value is in good agreement with theoretical estimates from 
percolation theory (a k, 0.686 in 2d (0)). 

A clear result from Fig. [T]A. is the absence of plateaus at crowding- 
dependent values for intermediate times and subthreshold obstacle densities. 
In fact, the only convincing plateau is obtained at the percolation threshold. 
This observation is confirmed by examination of the inset of Fig. [T]A Bot- 
tom, that presents aeff(i) for refined simulations (increased total simulation 
time, improved statistics) close to the threshold. These data show that one 
has to be very close to the threshold, i.e. the value of Q must be less than 
1% of its value at the threshold, and wait for a rather long time (i ^ 5 ms), 
in order to obtain a convincing power law regime. When the excluded frac- 
tion Q < 0.99^0 the curves fail to evidence plateaus, but instead show a 
concave upward curvature. We thus conclude that Fig. [T]A does not support 
the hypothesis of a power-law behavior with crowding-dependent exponents 
during the intermediate anomalous regime. 

To test the robustness of this conclusion, we ran simulations in dif- 
ferent conditions. Figure El] in the Supplementary Information section 
shows simulations in 2d but with various sizes of the tracer molecule (r^ = 
0.0 or 0.5 nm) and larger immobile obstacles (robs = 2.0 to 5.0 nm). The 
obtained results are very similar to that of Fig. [T]A. Similar results were also 
obtained when the random walks evolve on a square lattice (not shown). 
We also ran three-dimensional simulations (Fig. [Tf3). Because of the pro- 
hibitive computational cost of such 3d simulations, we could not reach the 
same statistical quality as in the 2d simulations. In particular, the maximal 
mean travelled distance was closer to the length of the spatial domain. De- 
spite these limitations, the conclusions are again similar to those brought 
about by the previous simulations. Our estimate for the anomalous exponent 
(a = 0.50(3)) is in good agreement with the theory (a ?a 0.526 in 3d (0)), 
and the estimated percolation threshold {9c = 0.87(1)) also agrees with the 
value estimated in Ref. (j36l ) (0.84). More importantly, these 3d simulations 
confirm the absence of local slope plateaus at intermediate times. 

Taken together, these results show that hindered diffusion by immo- 
bile obstacles does not feature transient anomalous regimes with crowding- 
dependent exponents. Instead, in agreement with theoretical predictions, 
obstacle density only modifies the crossover times between Brownian and 
transient anomalous regimes. Furthermore, for biologically realistic sizes, 
the simulations do not even show trace of real power law behavior at all 
during the transient anomalous regime. Indeed, the dynamics in the inter- 
mediate regime is completely dominated by the crossovers from and back to 
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Brownian diffusion. 



"Emulating" experimental results 

The extensive simulations reported above indicate that the intermediate 
anomalous power law expected near the percolation threshold may actually 
not be observable, in opposition to several experimental reports. A first 
reason for this discrepancy could be that the conditions of experimental 
measurements are far from the ideal situation of the above simulations. To 
test this hypothesis, we ran simulations of hindered diffusion with "much 
lower quality" than our previous ones, so as to test conditions closer to 
reported experiments (section "Emulating experimental results" in the Sup- 
plementary Information). The results from these simulations (Figure IS2p 
are indeed so imprecise that they can easily be misinterpreted in the sense 
of the existence of short pseudo-scaling laws. 

Even though, strictly speaking, one cannot compute a proper anomalous 
exponent in the transient regime (because no power law regimes are found), 
our simulation results can be used to estimate an approximate exponent 
during the anomalous regime, in order to, for instance, compare our results 
with those of works discussing how so-called anomalous scaling exponents 
vary. To this aim, we fitted the (oeff vs t) curves (cf Fig. [T]) by a nonlinear 
(log- normal) function to evaluate the minimal value of aeff(i) on the simu- 
lated time interval (see Methods). This value was then used as a "rational" 
estimate of the approximate exponent, a. 

First of all, the measurements of a on our data, in Fig. [2l shows that 
when plotted against the relative distance to the threshold e = 1 — 9 /6c (^c is 
the value of 6 at the percolation threshold), the data for 2d simulations and 
various molecule sizes collapse on a single behavior. This collapse suggests a 
generic behavior for approximate exponent that would only depend on space 
dimension. Moreover, the 2d and 3d data could both be fitted by the same 
nonlinear function (the lines in Fig. [2j): 



For 2d data, we obtained i^^^d) ^ o.75 ± 0.10 and af ^ = 0.69 ± 0.02 
and K^^'^^ = 0.23 ± 0.03 and a§'^^ = 0.49 ± 0.01 for the 3d case. Note that 
the value of a at the percolation threshold (e = 0) is given by the parameter 
ao and is expected to correspond to the theoretical values predicted by 
percolation theory. In this perspective, our estimates for ao are in very good 
agreement with theoretical values (a'^^'^^ ~ 0.686, a^^'^^ ~ 0.526 (2)). Hence, 
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we think that Eq. [2] yields an accurate functional form for the dependence 
of the approximate anomalous diffusion coefficient on crowding density (see 
the Discussion section for comparison with experimental results). 



Diffusion hindered by moving obstacles 

The previous section showed that poor statistics and finite size effects en- 
countered in typical experimental or numerical contexts may generate the 
appearance of pseudo-scaling laws with crowding-dependent exponents in 
the anomalous regime. This explains many of the disagreements between 
theory on the one hand and experimental and simulation reports on the 
other hand. However, several points remain unexplained, such as the re- 
ports of anomalous exponent values that are much smaller than the value 
expected at the percolation threshold (13) (see Introduction), although the 
latter is predicted to be the minimal accessible value. In this section, we 
examine the infiuence of obstacle mobility on the anomalous regime. 

Figure [3]A show simulations similar to those of Fig.[T]A with obstacle den- 
sity close to the percolation threshold (excluded volume fraction 6 = 0.418) 
but where obstacles move by Brownian diffusion. The reference curve in this 
figure is the red curve, that corresponds to .Dobs = i.e. immobile obstacles. 
The diffusion coefficient of the obstacles then progressively increases up to 
Dohs/D = 568, where D is the diffusion constant of the tracer molecules. 
The anomalous regime observed at long times with immobile obstacles (red 
curve) disappears as soon as Dohs/D > 1.0 (green curve). Hence, these re- 
sults show that anomalous diffusion is not expected to be observed at long 
times scales (around 0.1 ms) as soon as the obstacles diffuse, even with dif- 
fusion coefficient comparable to that of the tracer molecule. Note that this 



point was already suggested in Ref. (j3d ). and partly in Ref. (j37l ). 

However, diffusive obstacle movements are not the only possible ones. 
Figure [3f3 shows simulations with obstacles that fiuctuate around their equi- 
librium position. Note that in this figure, curves of the same color code for 
movements with identical mobility i.e. the same values of cr^^, (see Meth- 
ods). Comparing Fig. [3]A & B, it is clear that, with fluctuating movements, 
anomalous diffusion is more robust. In fact, with fluctuating obstacles, 
diffusion remains fully anomalous up to mobilities asr ~ 1.0 nm which cor- 
responds to Do\ys/D = 5.7 in the case of diffusive obstacles, i.e. a condition 
where the tracer movement is mostly Brownian. Thus, these results suggest 
that the transient diffusion regime could remain anomalous when obstacles 
fiuctuate around equilibrium positions, even when these fiuctuations are 
much larger than the mobility of the tracer molecule. 
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Another interesting property of anomalous diffusion with fluctuating ob- 
stacles is that when jitter amplitude is large enough, the fluctuations in effect 
remove the saturation of the MSD above the percolation threshold. Figure H] 
shows MSD plots similar to those in Fig. [3f3 except that the density of fluc- 
tuating obstacles is varied (while the fluctuation intensity is kept constant 
at a rather high value of asr = 1.0 nm). The top curve (the lavender one) 
corresponds to the same obstacle density as in Fig. [3B- Obstacle density is 
then progressively increased in the other curves of the figure, to values that 
would yield excluded fractions ranging from 6 = 0.418 to 9 = 0.762, if the 
obstacles were immobile. In other words, the density of obstacles used in 
Fig. |4]would be well above the percolation threshold with immobile obstacles 
{9c = 0.441, see Fig. [TJ^). At first inspection. Fig. U] shows that with very 
mobile fiuctuating obstacles, the anomalous diffusion regime is preserved for 
all curves, i.e. even when the obstacle density is much larger than the perco- 
lation threshold. However, the first crossover time — from the initial normal 
diffusion to the anomalous regime — increases with the obstacle density. 

Furthermore, a striking feature of the anomalous regimes obtained in 
suprathreshold conditions in Fig. H]is that the obtained values of the anoma- 
lous exponent a (i.e. the value of the local exponent a^s at the plateau) 
clearly decrease when the excluded volume fraction increases in the suprathresh- 
old domain. The values of these anomalous exponents in the suprathreshold 
domain are even found significantly below the value expected at the perco- 
lation threshold. With 6 = 0.709 for instance (the bottom curve of Fig. 2]), 
the value of a that would be deduced from the local slope for t > 1 ms is 
~ 0.40, which is much lower than theoretical value in 2d (~ 0.69). We thus 
propose that the presence of fluctuating obstacles in the suprathreshold do- 
main may explain some of the experimental reports of crowding-dependent 
anomalous exponents a as well as the fact that some of these values may be 
smaller than their theoretical value at the percolation threshold. 

Trapping time distributions 

One way to interpret hindered diffusion is to consider that the obstacles form 
cages that temporarily trap the tracer molecules. When the obstacles are 
mobile, they can be considered to be forming dynamic cages whose lifetimes 
are expected to be broadly distributed at high obstacle densities (jsl)- If this 
distribution is algebraic (and if the distribution exponent /3 G [1.0 — 2.0]), 
one generically expects anomalous diffusion, through a mechanism that is 
essentially the same as that of the CTRW model p, fl) (see Introduction) . 
In particular, in this case, the anomalous diffusion exponent a would be 
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expected to vary continuously with /? (q = /3 — 1). 

Because the traps formed by the moving obstacles are constantly reor- 
ganizing, it is difficult to define them. We opted for arguably the simplest 
definition. We considered that a walker is trapped by an obstacle as soon 
as the distance between them is less than some value rtrap (i-e. the distance 
between their centers is less than rtrap + ?'w + ^obs)- We thus measured the 
trapping time distribution of 2d tracer molecules in presence of robs = 5 nm 
obstacles that are either immobile (Fig. [5]A), fluctuating (Fig. [5p) or dif- 
fusive (Fig. IS3p . In each case, the first column presents conditions where 
the density of obstacles is low and diffusion remains Brownian throughout 
the simulation, whereas in the second column obstacle density and mobility 
were set so as to obtain a developed anomalous diffusion regime. 

As expected, with small densities of immobile obstacles (Fig. [5Kl) and 
for small values of the trapping distance rtrap, the trapping time distribu- 
tion is exponentially distributed, reflecting the correlation time of Brown- 
ian motion. For larger rtrap values, the distribution tends to become bi- 
exponential. When immobile obstacles are the source of the anomalous 
diffusion regime (Fig. [5K2), the distributions become broader and exhibit a 
power law tail with exponent around —3.5. With fluctuating obstacles, trap- 
ping time distributions with low obstacle densities (thus yielding Brownian 
diffusion. Fig. [5fBl) are very similar to those obtained with large obstacle 
densities (which yield anomalous diffusion. Fig. [5B2). In particular, in both 
case, the distributions display power-law tails with exponents around —3.8. 
Finally, in the case of diffusing obstacles, the distributions are qualitatively 
similar to the immobile case when diffusion is Brownian (Fig. IS3I A1). but 
with anomalous diffusion the power law-distributed long time tails observed 
above are much less obvious (Fig. [S3K21. 

Although we do not fully understand the various regimes observed, sev- 
eral interesting conclusions can already be drawn from these results. First 
of all, as expected, the distributions do broaden when obstacle density is 
large and diffusion anomalous. In certain cases, the distributions even show 
pronounced power-law tails. But we note that (i) in the case of mobile ob- 
stacles (either diffusing or fluctuating), the power law tails in the trapping 
time distributions are not specific to anomalous diffusion and (ii) the ob- 
served exponents /3 are found around 3.5 — 4.0, whereas in the CTRW model, 
anomalous diffusion is expected only for /3 G [1.0,2.0]. 

We thus conclude from these results that whatever the type of obsta- 
cle mobility, anomalous diffusion due to obstacles is not equivalent to a 
CTRW model, where anomalous diffusion would be explained by power- 
law distributed hopping times. Hence, this equivalence cannot be evoked 
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to explain the experimental observation of crowding-dependent anomalous 
diffusion exponents. 



Discussion 

The simulations reported here first of all confirm that hindered diffusion 
due to random immobile obstacles does not feature asymptotic anomalous 
scaling regimes where the exponent value would decay with obstacle density. 
The anomalous exponent is unique, only observable very near the percola- 
tion threshold, and obstacle density only modifies the crossover times be- 
tween Brownian and anomalous regimes. Moreover, for biologically realistic 
parameters, the MSD dynamics is completely dominated by the crossovers 
from and back to Brownian diffusion. 

One could first raise the issue of the realism of our simulation algorithms. 
We have used a simplified Monte-Carlo algorithm for molecule movement in 
a lattice-free space (see Methods) . For a subset of the conditions illustrated 
in the paper, we have compared our results with simulations obtained using 
a more classical Brownian dynamics algorithm (where the distance travelled 
at each time step is a Gaussian random variable) or with simulations of dif- 
fusion on a lattice. We have not found differences in the (log MSD vs log t) 
curves between these simulations and the ones illustrated above (not shown). 
More realistic diffusion models can be used, in particular using a real po- 
tential (e.g. Lennard-Jones) to simulate excluded volume or more realistic 
obstacle size distributions. Such an effort is accessible with modern day 
computing resources, albeit the total simulated real time is much smaller 
than in our paper (20 /is in. McGuffee and Elcock (jsil ) versus 10 — 100 ms in 
our case). Note however that, on their very short time scales, these detailed 
simulations are in general agreement with ours concerning the anomalous 
diffusion regimes (compare e.g. Figure 3B in McGuffee and Elcock (jssi ) with 
the bottom panels of Fig. [T]in the present work). Hence, we think that our 
results are generic and robust to refinements of the model. 

Yet, numerous experimental studies, in vitro as well as 



m VIVO , Esl ). have consistently reported anomalous diffusion regimes 
with variable values of the exponent a. A first explanation could reside 
in experimental limitations. Indeed, three major methods have been used 
to measure molecular diffusion in these conditions : single-particle tracking 
(SPT), fiuorescence recovery after photobleaching (FRAP) and fiuorescence 
correlation spectroscopy (FCS). Experimental errors in the determination of 
particle positions in SPT experiments has been shown to lead to spurious 
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(subdiffusive) power law regimes at short times (|38l). FRAP studies usually 
are limited in measurement time, which impedes the analysis of the long-tail 
phenomena expected in anomalous diffusion (figl ). FCS is also problematic 
as it rests on the time-correlation function of the fluorescent light emitted 
from the proteins, whose functional form is known for Brownian diffusion. 
The analytical formulation of the FCS correlation function is not available 
for anomalous diffusion, so that an adhoc function is usually employed, that 
lacks real theoretical support (^) (albeit recent advances could improve this 
caveat ([E))- Hence, even in in vitro and carefully controlled conditions, the 
experimental determination of some a value is very delicate. Moreover, inde- 
pendently from the measurement methods, the quantification of a in living 
cells could be influenced by intrinsic limitations due to the characteristics of 
the cell itself. In particular, our simulations (Fig. IS2p show that the finite 
size of the cells and/or the low statistics inherent to in vivo measurements, 
could by themselves give rise to spurious power law-like regimes, that would 
naturally be interpreted as crowding-dependent anomalous exponents. 

Another possible explanation could reside in the mobility of the obsta- 
cles. Indeed, our simulations confirm the conclusion of previous studies that 
the long-time anomalous regime typical of immobile obstacle disappears very 
rapidly as soon as the obstacles are mobile. However, our finding that diffu- 
sion amidst fiuctuating obstacles gives rise to extended anomalous regimes 
has interesting implications. First, the case of fluctuating obstacles seems 
rather relevant in the context of biological cells. Fluctuating obstacles could 
arise when they self-organize into noncrystalline but organized supramolecu- 
lar structures such as gel-like networks. This includes of course intracellular 
networks such as the actin or spectrin ones. For instance, a recent study 
measured the amplitudes of individual (i.e. not population averaged) MSDs 
of tracer particles within cytoplasms of eucaryotic cells (js^) and showed they 
need circa 10s to decorelate, thus suggesting intracellular rearrangements of 
the cytoskeleton or internal macromolecular complexes. This is thus an indi- 
rect support for the existence of a long range spatial self-ordering that could 
lead to fluctuating obstacle structures. But fiuctuating obstacles may be 
even more generically encountered. For instance, light-harvesting complexes 
of photosynthetic membranes are large-size obstacles that occupy between 
70 and 90% of the membrane area (jSOl ). Recent experimental observations 
revealed that their mobility in the membrane indeed consists of fluctuations 
around an equilibrium position, with fluctuation amplitudes that depend on 
the considered region in the membrane (j4ol ). 

Interestingly, our results with fluctuating obstacles could also explain 
some of the experimental reports obtained in vitro for tracers diffusing in 
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concentrated mobile crowders. Horton et al. ((51) for instance, studies 2d dif- 
fusion on artificial lipid bilayer membranes with avidin as both the crowder 
and the tracer. This study shows that even at small densities, avidin self- 
organizes into a spatially heterogeneous (possibly gel-like) structure, that 
persists (i.e. does not mix) for several tens of minutes. This type of spa- 
tial heterogeneity is typically expected to yield fluctuating obstacles. In 
agreement with our simulation results for fluctuating obstacles (Fig. [3|3), 
the authors find clear evidence for transient anomalous diffusion. Their 
measurements are also in striking agreement with our simulation results. 
For instance, the evolution with obstacle density of their estimates for the 
anomalous exponent (their Figure 5) has a shape very similar to our figure 
of the evolution of the approximate exponent a with the excluded fraction. 
Fig. [2] (which is not the case for the functional dependence proposed with 
mobile soluble crowders (f2oi)). Moreover, the estimation of the minimal 
value of the anomalous exponent in Horton et al. ^) (0.68) is also in very 
good agreement with our estimates in 2d (0.69). We consider that this very 
good correspondence between these experimental results and our simula- 
tions is a confirmation of the robustness of our prediction for fluctuating 
obstacles. 

Finally, we also presented evidence that whatever the type of mobility 
(or immobility), hindered diffusion due to obstacles is not equivalent to the 
CTRW model for trapping-based anomalous diffusion. The distribution of 
trapping-times due to the presence of obstacles broadens when the density of 
obstacles is large enough to yield anomalous diffusion, but it does not become 
a power-law with exponent in [1.0 — 2.0], which is a mandatory prerequi- 
site to yield anomalous diffusion in the CTRW model. Alternatively, the 
experimental observations of anomalous diffusion with crowding-dependent 
exponents could be explained by a real CTRW model, that is, it could arise 
if the movements of the tracer due to pure Brownian diffusion (i.e. in an ho- 
mogeneous space) were combined to weak reversible binding events to other 
molecular species. However, in this case again, the binding times must be 
power-law distributed with an exponent restricted to a narrow range. It is 
not clear how such a process could be obtained in a living cell. 

Hence, the major conclusion from this work is that, if one considers 
hindered diffusion due to obstacles as being the main source of anomalous 
diffusion within cells, the experimental observations of crowding-dependent 
anomalous diffusion are either spurious experimental results if obstacles are 
immobile or can be explained if the obstacles are present in high density 
and fluctuating around equilibrium positions within self-organized spatial 
structures. 
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Figure Legends 

Figure [H 

Evolution of the mean-squared displacement {Top Row) and cor- 
responding effective exponent Ocff {Bottom Row) during hindered 
diffusion by immobile obstacles in two (A) or three (B) dimen- 
sions. (A) 2d diffusion with excluded volume fraction ^ = 0, 0.06, 0.12, 
0.17, 0.22, 0.27, 0.31, 0.36, 0.40, 0.41, 0.42, 0.43, 0.45, 0.47 and 0.50 (from 
top to bottom). = 2.0 nm, robs = 5.0 nm, D = 4.4 fim^/s. Inset : Re- 
fined simulations with averages over 400 obstacle configurations and 9 = 
0.4345, 0.4365, 0.4380, 0.4400, 0.4405, 0.4410, 0.4415, 0.4430, 0.4450 (from 
top to bottom). Resulting estimates are 9c = 0.441(2) and a = 0.68(4). 
(B) 3D diffusion with 9 = 0, 0.567, 0.766, 0.843, 0.867, 0.874 and 0.876 
(from top to bottom), = 2.0 nm, robs = 10.0 nm, D = 7.7 /um^/s. These 
data yield 9c = 0.87(1) and a = 0.50(3). All data are averages over 10^ 
walkers per obstacle configurations and 10 obstacle configurations. Other 
parameters : 6t = 20 ns, 6r = 1.0 nm and = Ly = 5.0 ^m (2d) or 
Lx = Ly = Lz = 1-0 lira (3d). 

Figure [2], 

Variations of the approximate anomalous exponent a as a func- 
tion of excluded volume fraction 9. Space dimension, radius of the 
random walker and of the obstacles were {(i, rw,robs} = {2,0.0,5.0} (tri- 
angles), {2,0.5,2.0} (times signs), {2,2.0,5.0} (plus signs) or {3,2.0,10.0} 
(circles). The lines show nonlinear fits (see text for the function used and 
the parameter estimates) : the full (blue) line is a fit of all data points 
in 2d series as a whole, while the dashed (green) one is a fit of the 3d se- 
ries. The values of 9c were estimated from = f{t) curves (Fig. [TlandlSip. 

Figure [3l 

2D-diffusion hindered by mobile obstacles with various mobility. 

(A) Obstacles with diffusive motion. The obstacle-to-walker ratio of diffu- 
sion coefficients, Dohs/D = 0.0,5.7 x 10-^ 1.42, 5.68, 35.5, 142, 320 and 568 
(from bottom to top in the bottom panel). (B) Obstacles with fluctuating 
motion. The amplitude of the obstacle fluctuations are chosen identical to 
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that yielding diffusive motions in (A), i.e (from bottom to top in the bot- 
tom panel) asr = 0.0, 0.01, 0.50, 1.0, 2.5, 5.0, 7.5 and 10.0 nm. In both cases, 
the number of obstacles corresponded to 9 = 0.418 with immobile obstacles 
(i.e. close to the percolation threshold). Averages taken over 10^ walkers 
per obstacle configurations and 10 different obstacle configurations. Other 
parameters : D = 4.4 /im^/s, 5t = 20 ns and = Ly = 5.0 ^m. 

Figure [H 

Anomalous diffusion with fluctuating obstacles beyond the perco- 
lation threshold. Excluded volume fraction 9 = 0.418, 0.586, 0.646, 0.709 
and 0.762 from top to bottom, i.e. well above the percolation threshold 
with immobile obstacles (compare with figure [T]A). Obstacle mobility was 
asr = 1.0 nm. All other parameters were the same as figure [3j3. 

Figure [5l 

Trapping times distributions for hindered diffusion by immobile 
(A) or fluctuating (B) obstacles. In panels (Al, Bl), obstacle density 
is low enough [9 = 0.031) so that diffusion is Brownian throughout the 
simulation (see figure [Tl\) whereas in panels (A2, B2), obstacle density is 
close to the threshold {0 = 0.441) and diffusion is strongly anomalous. For 
mobile obstacles, mobility was chosen small enough to preserve anomalous 
behavior, i.e. asr = 1-0 nm for fiuctuating obstacles (B) (see figure [3]) . In 
each panel, the trapping distance rtrap = 0.5, 1.0, 2.50 and 5.0 nm (from 
top to bottom), respectively. Each curve is averaged over 10^ walkers per 
obstacle configurations and 100 obstacle configurations. The dashed lines 
locate power laws with exponent -3.5 (A2) or -3.8 (Bl-2). Other parameters: 
= 2.0 nm, Tobs = 5.0 nm, = 4.4 /zm^/s, 5t = 20 ns, = Ly = 
5.0 fim. 
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Supplementary Informations 



Emulating experimental results 

The extensive simulations reported above thus confirm the theoretical view- 
point of a unique anomalous exponent a whatever the intensity of molecular 
crowding. They also indicate that the expected intermediate anomalous 
power law may actually not be observable because it is often masked by the 
crossovers. Now, experimental data have consistently reported anomalous 
power laws at intermediate times with variable exponents (see Introduc- 
tion). How to explain such a discrepancy? A first possibility would be that 
the conditions of experimental measurements are far from the ideal situa- 
tion of the above simulations. To test this hypothesis, we ran simulations 
of hindered diffusion with "much lower quality" than our previous ones, so 
as to test conditions closer to reported experiments. In particular, we took 
into account the fact that the spatial domain in real experiments is usually 
limited by the size of the cell or the subcellular compartment in which the 
molecule diffusion is restricted. The quality of the statistics of the data in 
the experiments is also expected to be lower and usually, the sampling of 
the data points in time is much lower. 

Figure [S2] shows three-dimensional simulations in a space domain whose 
volume is one hundredth that of the one used in Fig. [TB above. In fact, the 
maximal mean travelled distance in Fig. [S2] is more than 5-fold the length 
of the space domain (we used periodic boundary conditions). We moreover 
averaged the data over a number of molecule per run that is ten times lower 
and a number of run that is only one third that of Fig. [TB- Finally, sam- 
pling frequency of the simulation to plot the results was also divided by ~ 50, 
yielding only 12 data points per curve. This figure illustrates the fact that 
the conjugation of poor statistics with finite size effects leads to the appear- 
ance of regimes in the local exponent curves that look very close to horizontal 
(Fig. IS21 Bottom). Importantly, the value of the local exponent (i.e. the 
value at the plateau) during these pseudo-scaling laws generically decreases 
with increasing excluded volume fraction, in agreement with several exper- 
imental reports. Note finally that these pseudo-scaling regimes usually last 
less than 1.5 decades and contain 3-5 data points. They are thus very sim- 
ilar to the pseudo-power laws with crowding-dependent exponents that are 
reported in several experimental as well as numerical studies. Hence these 
results show that the recurrent experimental reports of crowding-dependent 
anomalous diffusion exponents could be spurious effects due to a conjugation 
of finite size effects and poor statistics (w.r.t. simulations). 
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Even though our results show that the vs t relation does not 

evidence a power law regime (ie a plateau in the aeff(^) vs t curve), a tenta- 
tive approximate exponent, a, during the anomalous regime can be defined 
in order to, for instance, compare our results with those of works discussing 
how so-called anomalous scaling exponents vary. 

Our objective is to use the minimal value of the Oeff (t) vs t curve. Unfor- 
tunately, these curves are usually noisy, with significant "oscillations" super- 
imposed to the overall envelop (see e.g. Figure [H Bottom row), which pre- 
cludes naive attempts to use their minimal value as reliable estimates. But 
these curves (especially their concave part) can nicely be fitted by an inverted 

Log-Normal function of the form Qefr(t) = 1 — Aexp — {ln{t/tmin) /w)^ . 
The estimated value of A then yields the minimum value of the function on 
the interval, i.e. a = 1 — A. q is thus used here as a "rational" estimate of 
the approximate exponent. 

First of all, the measurements of a on our data, in Fig. O shows that 
when plotted against the relative distance to the threshold e = 1 — 9 /9c {9c is 
the value of 9 at the percolation threshold), the data for 2d simulations and 
various molecule sizes collapse on a single behavior. This collapse suggests a 
generic behavior for approximate exponent that would only depend on space 
dimension. Moreover, the 2d and 3d data could both be fitted by the same 
nonlinear function (the lines in Fig. [2]): 

d (e) = do + (1 - ao) (1 + i^) 



K + e 



For 2d data, we obtained i^^^d) ^ o.75 ± 0.10 and qJ,^'^^ = 0.69 ± 0.02 
and K'-^'^^ = 0.23 ± 0.03 and df,^''^ = 0.49 ± 0.01 for the 3d case. Note that 
the value of d at the percolation threshold (e = 0) is given by the parameter 
do and is expected to correspond to the theoretical values predicted by 
percolation theory. In this perspective, our estimates for ao are in very good 
agreement with theoretical values (a*-^'^^ ~ 0.686, a^^"^^ ~ 0.526 (0))- Hence, 
we think that Eq. [2] yields an accurate functional form for the dependence 
of the approximate anomalous diffusion coefficient on crowding density (see 
the Discussion section for comparison with experimental results). 
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Figure SI: Evolution of the mean-squared displacement (Top Row) 
and corresponding effective exponent Oeff (Bottom Row) during 
hindered diffusion by immobile obstacles in 2D. (A) 2d diffusion of 
small tracer molecules (r„ = 0.5 nm, robs = 2.0 nm) with excluded volume 
fraction 9 = 0, 0.10, 0.18, 0.26, 0.33, 0.39, 0.45, 0.48, 0.51, 0.52, 0.53, 0.58 
(from top to bottom). Corresponding estimates are 6c = 0.52(2) and a = 
0.64(5). (B) 2d diffusion of point-like tracers (r„ = 0.0 nm, robs = 5.0 nm) 
with (9 = 0, 0.27, 0.47, 0.61, 0.62, 0.63, 0.648, 0.654, 0.545, 0.72 and 0.79 
(from top to bottom); yielding 6c = 0.66(1) and a = 0.69(3). All data are 
averages over 10^ walkers per obstacle configurations and 10 (A) or 100 (B) 
obstacle configurations. Other parameters : D = 4.4 /im^/s, St = 10 ns (A) 
or 7.7 fiTo^/s and 20 ns, respectively (B), 6r = 1.0 nm and = Ly = 5.0 /um 
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Figure S2: Simulations of hindered diffusion by immobile obstacles 
in poor conditions (large finite size eff"ect, low statistics, poor sam- 
pling) . The simulation conditions are identical to figure [TjB, i-e. d = 3, 
= 2.0 nm, Tobs = 10.0 nm except = Ly = = 0.10 ^um, averaging 
over 10^ walkers per obstacle configuration and 3 obstacle configurations 
and circa 50- fold less data points than in figure [TJ3 (open squares). = 0, 
0.448, 0.630, 0.773, and 0.851 (from top to bottom). Lines are guides to the 
eyes. 
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Figure S3: Trapping times distributions for hindered diffusion by 
immobile diffusive (B) obstacles. Same simulation parameters as in 
Figure [5] except that the obstacles move through Brownian diffusion with 
diffusion coefficient Dq-^^/D = 5.7 x 10~^. All other parameters as in Fig- 
ure O 
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